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Calculation of current and order parameter distribution in inhomogeneous superconduc- 
tors is often based on a self-consistent solution of Eilenberger equations for quasiclassical 
Green's functions. Compared to the original Gorkov equations, the problem is much sim- 
plified due to the fact that the values of Green's functions at a given point are connected 
to the bulk ones at infinity (boundary values) by "dragging" along the classical trajecto- 
ries of quasiparticles. In finite size systems, where classical trajectories undergo multiple 
refiections from surfaces and interfaces, the usefulness of the approach is no longer ob- 
vious, since there is no simple criterion to determine what boundary value a trajectory 
corresponds to, and whether it reaches infinity at all. Here, we demonstrate the modifi- 
cation of the approach based on the Schophol -Maki transformation, which provides the 
basis for stable numerical calculations in 2D. We apply it to two examples; generation 
of spontaneous currents and magnetic moments in isolated islands of d-wave supercon- 
ductor with subdominant order-parameters s and dxy, and in a grain boundary junction 
between two arbitrarily oriented d-wave superconductors. Both examples are relevant to 
the discussion of time-reversal symmetry breaking in unconventional superconductors, 
as well as for application in quantum computing. 



1 Introduction 

Pairing symmetry of unconventional superconductors can produce time-reversal 
symmetry breaking states. El Especially interesting are higJi Tc cuprates, with their d- 
wave symmetry, since the recently developed technologyB allows fabrication of such 
structures, with controllable characteristics, as 7r-junctions, submicron size "(/)o"- 
junctionsH (i.e. junctions with equilibrium phase difference (f>o which is neither Q 
nor 7r),Q tt-SQUIDsB or 7r/2-SQUIDs, □ and superconducting qubit prototypes. Ll 
Therefore, quantitative prediction of properties of such restricted systems becomes 
relevant not just from the theoretical point of view. 

General approach to such calculations is based on Gorkov equations for Green's 
functions of the superconductor. It is limited only by the applicability of BCS-like 
"mean field" description of the superconducting state. This description is now 
considered valid on the phenomenological level, B independently from further devel- 
opments in the first principles' theory of high Tj, superconductivity. Quasiclassical 
limit of these equations, Eilenberger equations, u is strictly speaking valid only when 
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|A| ^ Ep,. Unlike in conventional superconductors, this condition is not completely 
satisfied in high Tc superconductors but is still a good approximation. 



2 Quasiclassical Approach 

We use the standard approach based on Eilenberger equations i for quasiclassical 
Green's functions 

VF • V5 + [u;?3 + A,g]= 0, (1) 
with normalization condition 7j ^ ~ 1, where lu is the Matsubara frequency and 

The matrix Green's function g and the superconducting order parameter A are both 
functions of the Fermi velocity vp and position r. A is determined by the (2D) 
self-consistency equation 



A{e)^27:N{0)Tj2{yee'f{0'))e 



LJ>0 

where 9 is the angle between -vp and the a;-axis, Vggi interaction potential, A'^(O) 
density of states at the Fermi surface, and {■■■)g represents averaging over 9. Gen- 
erally, it is possible to obtain a mixture of different symmetries of the order param- 
eter, e.g. A = Ax2_y2 + Axy + As where Ax2_y2 = Aicos20, A^y = A2sin2^, 
and As are the dominant (1^2 _y2 component, and the subdominant d^y and the 
s components of the order parameter respectively. The corresponding interaction 
potential, Veo' — Vdi cos 20 cos 2^^' -I- Vd2 sin 20 sin 20' -I- Vs, must be substituted in 
the self-consistency equation. The current density j (r) is found from g as 

j^-47:ieN{0)Tj2{^Fg)g 

For numerical calculations, it is conventional@00 to parameterize the quasi- 
classical Green's functions by so called coherent functions a, b via 

I — ab 2a , . 

' = TT^b^ ■^ = TT^b- 

Functions a and b satisfy two independent, but nonlinear, equations 

Vp ■ Va = A — A*a^ — 2a;a 
-wp . V6 A* - Ab^ - 2ujb. (3) 

From these equations it follows that a(~vp) — b*{vp) and b{—vp) = a*{vp). One 
should solve these equations along all possible quasiclassical trajectories and perform 
the summation over the trajectories to calculate the order parameter current density. 
To find a and b along the trajectories, one needs to use boundary conditions at the 
ends of the trajectories. In macroscopically large systems, one usually assumes that 
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Figure 1: Absolute values of the dominant (left) and subdominant (right) order parameters for 
a small square d-wave super conductor in the presence of subdominant s-wave order parameter. 
The orientation of the main order parameter is 45° rotated with respect to the boundaries. 



all the trajectories go deep into the bulk of the superconductor, i.e., it is possible 
to use the bulk solutions 



A 



a± 



A* 



(4) 



with Q — ^/uJ^~{~\A\^ , as the boundary conditions at infinity. In the case of "re- 
stricted" systems, the above assumption is no longer self-evident. Nevertheless, we 
will sec that a stable numerical procedure can be still developed. 

Numerical calculation is stable if the integration for a (b) is taken in the (op- 
posite) direction of v^?. When A is a constant, the solution of Eq. (0) for a is 



1 + —{a-i - a 
il I a j — a 
A 



a* 



for fir > 1 



(5) 



where and af are the values of a at the initial (r^) and final (r/) points of the 
trajectory, and t = \rf — ri\/vF is proportional to the distance between the initial 
and final points along the trajectory. It is clear that the solution for a relaxes to 
the bulk value a+ at the distance L = \-p/2il which is of the order of the coherence 
length ^g. In other words, when the quasiparticle moves away from the initial point 
at a distance of a few ^o's, any information about the initial point becomes lost. 
This observation is also valid for b, and is crucial for what that follows. 

Let us now specifically consider a restricted system. After integrating over a few 
^o's, af will be almost independent of Oj, although it may never coincide with the 
bulk value a+ . This solution corresponds to a simple exponential relaxation of the 
functions a and b to their local "steady-state" values defined by the local value of 
the order parameter. This value is the limit for the functions a and 6 at this spatial 
point. Such relaxation of Green's functions significantly simplifies the numerical 
solution of the self-consistent two-dimensional problem. The system therefore has 
no memory of the local values of A beyond several along the trajectory. 

In order to calculate a, we define all possible trajectories going through a given 
point of the system. Along each of them we move back (in the direction opposite 
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to vp) a cutoff distance (about f 0^o~ 20^o) and choose that point as the beginning 
of the trajectory. We set the bulk solution (0+) as the initial value for a at that 
point (this really does not matter too much because the system has no memory) and 
integrate along the trajectory, taking into account the reflections at the boundaries, 
until we get back to the calculation point. Calculation for b is the same, except that 
the direction of integration is now opposite to v^. These calculations are repeated 
for all trajectories and all points of the system (on a mesh). After each iteration, 
we calculate the new A and use it for the next iteration until the self-consistency 
is achieved. We found this method to be very stable and independent of the value 
of the "cutting" distance. One should note that the above procedure is not vahd in 
the presence of magnetic field, because a path dependent phase will be accumulated 
to the Green's functions, and the above mentioned relaxation mechanism along the 
trajectory does not work anymore. 

3 Results 

To illustrate the approach, we performed self-consistent calculations of the order 
parameter in a small (20^o x 20^o) square region of c?-wave superconductor in the 
presence of a subdominant s-wave order parameter. The crystalographic a and b 
directions of the dominant order parameter make 45° angle with respect to the 
boundaries of the square. We used a random subdominant order parameter at the 
first iteration in order to avoid imposing any assumption on the phase of the second 
order parameter. 

The two components of the order parameter are displayed in Fig. 0. The dom- 
inant order parameter is suppressed at the boundaries of the square. This is due to 
the special orientation of the order parameter which requires all quasiparticles to 
face opposite sign of the order parameter after reflecting from the boundaries. □ As 
a result, the subdominant order parameter (s) with a tt/2 phase shift with respect 
to the dominant order parameter, is the main contributor at the boundaries. The 
spontaneous current distribution is displayed in Fig. |^. The current does not flow 
in the same direction at all the four edges, but changes the direction from one edge 
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Figure 3: Absolute value of the order parameter for a grain boundary junction between two d-wave 
superconductors. The grain boundary is a vertical line located in the middle. The orientation of 
the order parameter is 0° on the left and 45° on the right. 



to another closing its path towards the center of the square. The magnetic field 
produced by the current is also shown in Fig. |2| The maximum value of the mag- 
netic field is of the order of 10~^G. Notice that the direction of the filed changes 
from one edge to another. Thus, the total flux produced by such a field is zero. 

We also have calculated the spontaneous current and magnetic field distribu- 
tions in a d-wave grain boundary junction. The system consists of a finite square 
(30^0 X 30^0 ) rnade of d-wave superconductor divided into two equal parts separated 
by a grain boundary junction. The order parameter has dx2_y2 symmetry with 0° 
orientation on the left side of the grain boundary and 45° on the right side. We 
include a dxy subdomiiiant order parameter by adding an attraction potential in 
that channel (see Ref. □). The magnitude of potential is chosen in such a way ta 
have a transition temperature Tc2 = O.lTc (in the absence of the dominant order). □ 
We choose a phase difference of Ac/) = 7r/2 between the two sides. This actually 
corresponds to the equilibrium phase differexice of the junction at which the total 
current passing through the junction is zero.Q Calculations are done at T = O.OSTc. 

The results of calculation of the spontaneous current and magnetic field dis- 
tributions are displayed in Fig. |^. Notice that the current is not symmetric with 
respect to the grain boundary. On the left side (with 0° orientation), the current 
returns along the diagonal, whereas on the right side (45° orientation) it forms two 
vortices and antivortices near the edges— These vortices are a consequence of the 
chiral nature of the d -t- id' symmetry. □'Ej Magnetic field is peaked at the location 
of vortices with a maximum of the order of 10~^G. ft is important to emphasize 
that, unlike in the previous case, the existence of the spontaneous current in this 
system does not depend on the presence of a subdominant order parameter (al- 
though we assumed a subdominant component here). Addition of a subdominant 
order paramfiter will actually suppress the spontaneous current at the boundary 
(see Refs.yy). 
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4 Conclusions 



We described a method to calculate equilibrium properties in finite size supercon- 
ducting systems. We presented the results of our calculations for the distribution 
of spontaneous current and magnetic field in two systems: a small square region 
of a d-wave superconductor with a pair breaking boundary, and a d-wave grain 
boundary junctions between two differently oriented d-wave superconductors. The 
method described here is quite general and can be applied to any 2D geometry with 
proper boundary conditions. Presence of external magnetic field invalidates the 
method described here. The self-generated magnetic field due to the spontaneous 
currents, however, is usually very small so that its effect can be neglected in most 
calculations. 
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